/* ****************************************************************************
 *
 * Copyright (c) Microsoft Corporation. 
 *
 * This source code is subject to terms and conditions of the Microsoft Public License. A 
 * copy of the license can be found in the License.html file at the root of this distribution. If 
 * you cannot locate the  Microsoft Public License, please send an email to 
 * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound 
 * by the terms of the Microsoft Public License.
 *
 * You must not remove this notice, or any other, from this software.
 *
 * THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
 * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
 * PARTICULAR PURPOSE.
 * ***************************************************************************/


#if !defined(AFX_HISTOGRAM_H__DF74A789_2878_4CB0_BAA9_125272F0E7E1__INCLUDED_)
#define AFX_HISTOGRAM_H__DF74A789_2878_4CB0_BAA9_125272F0E7E1__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "VectorF.h"
#include <algorithm>

/*! \brief A histogram class

  This class computes the histogram of a vector.

\par Template parameters
	
	- T type of input data (can be any: float, double, int, size_t, etc...)
	- TOut type of output data: float or double. (default is double)
	- VectorT container (vector-like) of T
	- VectorTOut container (vector-like) of TOut

\par Definition
  
Let \f$\mathcal N\f$ a data set of size \f$n\f$ defined on a domain \f$[a,b]\f$, 
the histogram of \f$\mathcal N\f$ is the plot of the density function \f$\phi(x)\f$ of
\f$\mathcal N\f$. \f$\phi\f$ is such that
\f[ \int_{-\infty}^{+\infty}\phi(x)\mathrm dx = n \f]
	
The normalized histogram is the plot of the normalized density function \f$\tilde \phi(x)\f$ defined by
\f[ \tilde \phi(x) = \frac{\phi(x)}{\int_{-\infty}^{+\infty}\phi(x)\mathrm dx}\f]
so that
\f[ \int_{-\infty}^{+\infty}\tilde \phi(x) \mathrm dx = 1 \ . \f]
	  
In the case of a discretized density function (implemented here) and for a step \f$dx\f$,
the density \f$\phi_i\f$ in the i-th region is defined by
\f[ \phi_i = \frac{d_i}{dx} \f] 
where is the number of data entry \f$d_i\f$ such that
\f[ a+i \; dx \le d_i < a+(i+1)dx \f]
		
Hence, the integral of the discretized density function can be computed as
\f[ \int_{-\infty}^{+\infty}\phi(x)\mathrm dx = \sum_i \phi_idx = \sum_i d_i = n \f]
and the discretized normalized density function is
\f[ \tilde \phi_i = \frac{\phi_i}{n}\f]
so that
\f[ \int_{-\infty}^{+\infty}\tilde \phi(x)\mathrm dx = \sum_i \tilde \phi_idx = 1 \f]
		
\par Why do normalized histograms show high values for low number of data ?
		  
This is easily understood by considering the following:
\f[ \tilde \phi_i = \frac{d_i}{n dx} \f]

 If \f$n=1\f$ and \f$d_i=1\f$, \f$\tilde \phi_i = \frac{1}{dx}\f$. If \f$dx\f$ is lower than 1, the normalized histogram will
have high values !

\par Moments:

	The moments (average, standard deviation, skewness, etc.) are computed using 
the algorithm of the Numerical recipies (see Numerical recipies in C, Chapter 14.1, pg 613).

\par Example:

	This example shows the typical use of the class:
\code
	// a vector containing the data
	vector<float> vData;
	// Creating histogram using float data and with 101 containers,
	THistogram<float> histo(101);
	// computing the histogram
	histo.Compute(vData);
\endcode

Once this is done, you can get a vector with the histogram or the normalized histogram (such that it's area is 1):
\code
	// getting normalized histogram
	vector<float> v=histo.GetNormalizedHistogram();
\endcode

\par Reference
	
	Equally spaced acsissa function integration (used in #GetArea): Numerical Recipies in C, Chapter 4.1, pg 130.

  \ingroup AlgoToolsGroup
*/

template<
	typename T = float, 
	typename TOut = double, 
	class VectorT = std::valarray<T>,
	class VectorTOut = std::valarray<TOut> >
class THistogram  
{
public:
	//! \name Constructors
	//@{
	/*! Default constructor
	\param nCounters the number of histogram containers (default value is 10)
	*/
	THistogram(size_t nCounters = 10);
	virtual ~THistogram()				{ Clear();};
	//@}

	//! \name Histogram computation, update
	//@{
	/*! \brief Computes histogram of vector v
	\param v vector to compute histogram
	\param bComputeMinMax set to true if min/max of v have to be used to get the histogram limits	
	
	This function computes the histogram of v and stores it internally.
	\sa Update, GetHistogram
	*/
	void Compute( const VectorT& v, bool bComputeMinMax = true);

	/*! \brief Computes histogram of vector v
	\param v vector to compute histogram
	\param bComputeMinMax set to true if min/max of v have to be used to get the histogram limits	
	
	This function computes the histogram of v and stores it internally.
	\sa Update, GetHistogram
	*/
	void Compute( const gfxinterface::CVectorF& v, bool bComputeMinMax = true);
	//! Update histogram with the vector v
	void Update( const VectorT& v);
	//! Update histogram with t
	void Update( const T& t);
	//@}

	//! \name Resetting functions
	//@{
	//! Resize the histogram. Warning this function clear the histogram.
	void Resize( size_t nCounters );
	//! Clears the histogram
	void Clear()											{  m_vCounters.free();};
	//! resets counters to 0
	void ResetCounters();
	//@}

	//! \name Setters
	//@{
	/*! This function sets the minimum of the histogram spectrum. 
	The spectrum is not recomputed, use it with care
	*/
	void SetMinSpectrum( const T& tMin )					{	m_tMin = tMin; ComputeStep();};

	/*! This function sets the minimum of the histogram spectrum. 
	The spectrum is not recomputed, use it with care
	*/
	void SetMaxSpectrum( const T& tMax )					{	m_tMax = tMax; ComputeStep();};
	//@}

	//! \name Getters
	//@{
	//! return minimum of histogram spectrum
	const T& GetMinSpectrum() const							{	return m_tMin;};
	//! return maximum of histogram spectrum
	const T& GetMaxSpectrum() const							{	return m_tMax;};
	//! \brief returns the maximum of occurence
	long GetMaxOccurence() const;
	//! \brief returns the maximum of occurence
	long GetMinOccurence() const;
	//! \brief returns the maximum of occurence, normalized with respect to the area
	TOut GetNormMaxOccurence() const						{	return static_cast<TOut>(GetMaxOccurence()/static_cast<TOut>(GetArea()*m_dStep));};
	//! \brief returns the minimum of occurence
	TOut GetNormMinOccurence() const						{	return static_cast<TOut>(GetMinOccurence()/static_cast<TOut>(GetArea()*m_dStep));};

	//! return step size of histogram containers
	TOut GetStep() const									{	return m_dStep;};
	//! return number of points in histogram
	long GetSum() const										{	return m_lArea;};
	/*! \brief returns area under the histogram 

		The Simpson rule is used to integrate the histogram.
	*/
	long GetArea() const									{	return m_lArea;};

	/*! \brief Computes the moments of the histogram

	\param vData dataset
	\param fAve mean
	\f[ \bar x = \frac{1}{N} \sum_{j=1}^N x_j\f]
	\param fAdev mean absolute deviation
	\f[ Adev(X) = \frac{1}{N} \sum_{j=1}^N | x_j - \bar x |\f]
	\param fVar average deviation:
	\f[ \mbox{Var}(X) = \frac{1}{N-1} \sum_{j=1}^N (x_j - \bar x)^2\f]
	\param fSdev standard deviation:
	\f[ \sigma(X) = \sqrt{Var(\bar x) }\f]
	\param fSkew skewness
	\f[ \mbox{Skew}(X) = \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^3\f]
	\param fKurt kurtosis
	\f[ \mbox{Kurt}(X) = \left\{ \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^4 \right\} - 3\f]

	*/
	static void GetMoments(const VectorT& vData, TOut& fAve, TOut& fAdev, TOut& fSdev, TOut& fVar, TOut& fSkew, TOut& fKurt);

	//! return number of containers
	size_t GetSize() const											{	return m_vCounters.size();};
	//! returns i-th counter
	size_t operator [] (size_t i) const								{	ASSERT( i < m_vCounters.size() ); return m_vCounters[i];};
	//! return the computed histogram
	const std::valarray<long>& GetHistogram() const					{	return m_vCounters;};
	//! return the computed histogram, in TOuts
	void GetHistogramD(VectorTOut& v) const;
	/*! return the normalized computed histogram 

	\return the histogram such that the area is equal to 1
	*/
	void GetNormalizedHistogram(VectorTOut& v) const;
	//! returns left containers position
	void GetLeftContainers(VectorTOut& v) const;
	/*! \brief return origin, step and number of steps for left containers

	*/
	void GetLeftContainers( T& origin, TOut& step, size_t& nSteps)			{	origin = m_tMin; step = m_dStep; nSteps=m_vCounters.size();};
	//! returns center containers position
	void GetCenterContainers(VectorTOut& v) const;
	//@}
protected:
	//! Computes the step
	void ComputeStep()												{	m_dStep = static_cast<TOut>(m_tMax-m_tMin) / m_vCounters.size();};
	//! Data accumulators
	std::valarray<long> m_vCounters;
	//! minimum of dataset
	T m_tMin;
	//! maximum of dataset
	T m_tMax;
	//! Area
	long m_lArea;
	//! width of container
	TOut m_dStep;
};


template<typename T, typename TOut, class VectorT, class VectorTOut>
THistogram<T,TOut,VectorT,VectorTOut>::THistogram(size_t nCounters)
: m_tMin(0), m_tMax(0), m_dStep(0)
{
	m_vCounters.resize(nCounters,0);
}

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::Resize( size_t nCounters )
{
	Clear();

	m_vCounters.resize(nCounters,0);

	ComputeStep();
}

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::ResetCounters()
{
	for (size_t i=0;i<m_vCounters.size();++i)
		m_vCounters[i]=0;
	m_lArea=0;
}

template<typename T, typename TOut, class VectorT,class VectorTOut>
long THistogram<T,TOut,VectorT,VectorTOut>::GetMaxOccurence() const
{
	long m=0;
	for (size_t i = 0; i<m_vCounters.size();++i)
		m=__max( m, m_vCounters[i] );

	return m;
//	return std::max<long>( m_vCounters ); 
};

template<typename T, typename TOut, class VectorT,class VectorTOut>
long THistogram<T,TOut,VectorT,VectorTOut>::GetMinOccurence() const
{	
	if (m_vCounters.size()==0)
		return 0;

	long m=m_vCounters[0];
	for (size_t i = 1; i<m_vCounters.size();++i)
		m=__min( m, m_vCounters[i] );

	return m;
//	return std::min<long>( m_vCounters ); 
};

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::Compute( const VectorT& v, bool bComputeMinMax)
{
	using namespace std;
	size_t i;
	int index;

	m_lArea=0;
	if (m_vCounters.size()==0)
		return;

	if (bComputeMinMax)
	{
		m_tMax = m_tMin = v[0];
		for (size_t i =1;i<v.size();++i)
		{
			m_tMax = __max( m_tMax, v[i]);
			m_tMin = __min( m_tMin, v[i]);
		}

//		m_tMax = max<T>( v );
//		m_tMin = min<T>( v );

		if (algotools::IsEqual(m_tMax,m_tMin))
		{
			m_tMax=m_tMin+1;
		}
	}

	ComputeStep();

	for (i = 0;i < v.size() ; ++i)
	{
		index=static_cast<int>( floor( ((TOut)(v[i]-m_tMin))/m_dStep )) ;

		if (index >= static_cast<int>(m_vCounters.size()) || index < 0)
			continue;

		m_vCounters[index]++;
	}

	m_lArea=m_vCounters.sum();
}

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::Compute( const gfxinterface::CVectorF& v, bool bComputeMinMax)
{
	using namespace std;
	size_t i;
	int index;

	m_lArea=0;
	if (m_vCounters.size() == 0)
		return;

	if (bComputeMinMax)
	{
		m_tMax = m_tMin = v[0];
		for (size_t i =1;i<v.size();++i)
		{
			m_tMax = __max( m_tMax, v[i]);
			m_tMin = __min( m_tMin, v[i]);
		}

		if (algotools::IsEqual(m_tMax,m_tMin))
		{
			m_tMin=0;
			m_tMax=1;
		}

	}

	ComputeStep();

	for (i = 0;i < v.size() ; ++i)
	{
		index=static_cast<int>(floor( ((TOut)(v[i]-m_tMin))/m_dStep )) ;

		if (index >= static_cast<int>(m_vCounters.size()) || index < 0)
			continue;

		m_vCounters[index]++;
	}

	m_lArea=m_vCounters.sum();	
}

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::Update( const VectorT& v)
{
	if (m_vCounters.size()==0)
		return;

	ComputeStep();

	size_t uSize = m_vCounters.size();

	int index;
	for (size_t i = 0;i < uSize ; ++i)
	{
		index = static_cast<int>(floor(static_cast<TOut>(v[i]-m_tMin)/m_dStep));

		if (index >= static_cast<int>(m_vCounters.size()) || index < 0)
			continue;

		m_vCounters[index]++;
	}
	m_lArea=m_vCounters.sum();	
}

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::Update( const T& t)
{	
	int index=static_cast<int>(floor( static_cast<TOut>((t-m_tMin)/m_dStep) ) ) ;

	if (index >= static_cast<int>(m_vCounters.size()) || index < 0)
		return;

	m_vCounters[index]++;
	m_lArea=m_vCounters.sum();	
};

template<typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::GetHistogramD(VectorTOut& v) const
{
	v.resize(m_vCounters.size());
	for (size_t i = 0;i<m_vCounters.size(); ++i)
		v[i]=static_cast<TOut>(m_vCounters[i]/m_dStep);
}

template <typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::GetLeftContainers(VectorTOut& v) const
{
	v.resize( m_vCounters.size());

	for (size_t i = 0;i<m_vCounters.size(); ++i)
		v[i]= m_tMin + i*m_dStep;
}

template <typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::GetCenterContainers(VectorTOut& v) const
{
	v.resize( m_vCounters.size());

	for (size_t i = 0;i<m_vCounters.size(); ++i)
		v[i]= static_cast<TOut>(m_tMin + (i+0.5f)*m_dStep);
}

template <typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::GetNormalizedHistogram(VectorTOut& v) const
{
	v.resize( m_vCounters.size());
	TOut dArea = (TOut)GetArea()*m_dStep;

	for (size_t i = 0;i<m_vCounters.size(); ++i)
	{
		v[i]= static_cast<TOut>(m_vCounters[i]/dArea);
	}
};

template <typename T, typename TOut, class VectorT,class VectorTOut>
void THistogram<T,TOut,VectorT,VectorTOut>::GetMoments(const VectorT& vData, TOut& fAve, TOut& fAdev, TOut& fSdev, TOut& fVar, TOut& fSkew, TOut& fKurt)
{
	size_t j;
	double ep=0.0,s,p;
	const size_t n = vData.size();

	if (n <= 1)
		// nrerror("n must be at least 2 in moment");
		return;

	s=0.0; // First pass to get the mean.
	for (j=0;j<n;j++)
		s += vData[j];
	
	fAve=(TOut)(s/(n));
	fAdev=fVar=fSkew=fKurt=0.0; 
	/* Second pass to get the first (absolute), second,
	third, and fourth moments of the
	deviation from the mean. */

		for (j=0;j<n;j++) 
		{
			fAdev += (TOut)fabs(s=vData[j]-(fAve));
			ep += s;
			fVar += (TOut)(p=s*s);
			fSkew += (TOut)(p *= s);
			fKurt += (TOut)(p *= s);
		}


	fAdev /= n;
	fVar=(TOut)((fVar-ep*ep/n)/(n-1)); // Corrected two-pass formula.
	fSdev=(TOut)(sqrt(fVar)); // Put the pieces together according to the conventional definitions. 
	if (fVar) 
	{
		fSkew /= (n*(fVar)*(fSdev));
		fKurt=(TOut)((fKurt)/(n*(fVar)*(fVar))-3.0);
	} 
	else
		//nrerror("No skew/kurtosis when variance = 0 (in moment)");
		return;
};

#endif // !defined(AFX_HISTOGRAM_H__DF74A789_2878_4CB0_BAA9_125272F0E7E1__INCLUDED_)
